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PREFACE 


This  study  is  an  attempt  to  apply  discrete  optimal  con- 
trol theory  to  the  design  of  a digital  flight  control  system. 

I have  chosen  the  terrain  following  task  to  illustrate  the 
feasibility  of  such  a control  system.  It  would  be  helpful 
to  the  reader  if  he  has  a basic  knowledge  of  optimal  control 
theory  and  aircraft  dynamics.  I believe  the  develpoment  of 
the  study  is  straight  forward  and  could  similarly  be  applied 
to  any  flight  control  task.  I feel  that  digital  control  will 
bring  about  a minor  revolution  in  the  realm  of  aircraft  per- 
formance. With  the  continued  improvement  in  reliability  and 
miniaturization  of  solid  state  devices,  the  flight  control 
system  of  the  future  will  be  easily  modified  through  soft- 
ware to  perform  a large  variety  of  tasks.  I further  feel 
that  discrete  optimal  control  theory  will  form  the  foundation 
for  the  development  of  the  required  control  laws.  Optimal  con- 
trol has  an  intrinsic  beaty  that,  unfortunately,  I have  only 
glimpsed. 

I would  like  to  express  my  deepest  gratitude  to  my  AFIT 
thesis  advisor.  Captain  James  E.  Negro,  for  his  encouragement 
and  many  helpful  suggestions.  I am  also  indebted  to  Michael 
J.  Breza  of  the  Aerospace  Systems  Division,  Deputy  for  Devel- 
opment Planning,  who  sponsored  this  research  effort.  Finally, 

I wish  to  thank  Major  James  E.  Funk  for  his  previous  work 
which  provided  me  with  an  optimal  path  for  the  terrain  fol- 
lowing mission. 

Ross  L.  Simmons 
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Abstract 

Through  the  use  of  state  space,  continuous  optimal  con- 
trol, and  discrete  optimal  control,  a digital  flight  control 
system  was  designed  for  the  terrain  following  task.  After 
formulating  the  aircraft  linear  perturbation  model,  the 
deterministic  regulator  problem  was  solved  with  a quadratic 
performance  index  to  provide  the  desired  continuous  closed 
loop  system.  The  system  and  performance  index  were  then  dis- 
cretized to  form  a discrete  deterministic  regulator  problem. 
This  discrete  regulator  problem  was  then  solved  as  a function 
of  sample  rate  using  eigenvector  decomposition  to  determine 
a minimum  acceptable  rate  for  sampling.  The  effects  of  sample 
rate  on  the  system  were  then  examined.  A sample  rate  of  five 
hertz  was  shown  to  be  high  enough  to  adequately  form  the  de- 
sired controls.  A reference  command  generator  based  on  con- 
stant energy  path  legs  was  developed  to  provide  the  required 
reference  states  and  control  inputs.  The  reference  terrain 
following  path  was  generated  by  an  optimal  cubic  spline  al- 
gorithm. The  aircraft  was  shown  to  track  the  desired  path 
in  a highly  acceptable  manner  through  the  use  of  a hybrid 
simulation.  The  design  method  utilized  is  recommended  for 
consideration  in  designing  the  digital  flight  control  system 
for  other  flight  control  tasks,  v 
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FOR  THE  ADVANCED  TACTICAL  FIGHTER 
USING  DISCRETE  OPTIMAL  CONTROL 

I . Introduction 

The  breakthroughs  of  the  past  decade  in  the  miniatur- 
ization and  improved  reliability  of  electronic  devices  have 
provided  the  control  engineer  with  an  alternate  means  of 
implementing  an  aircraft  flight  control  system.  Instead  of 
the  bellcranks,  pushrods,and  other  mechanical  linkages  nor- 
mally associated  with  an  aircraft  flight  control  system,  it 
is  now  possible  to  use  electrical  signal  paths  for  informa- 
tion transfer.  The  replacement  provides  a significant  de- 
crease in  the  weight  and  required  size  of  the  aircraft.  This 
implementation  has  been  accomplished  in  the  F-16  aircraft  by 
the  use  of  an  analog  fly  by  wire  flight  control  system.  Fur- 
ther breakthroughs  in  the  miniaturization  and  reliability  of 
electronic  devices  have  made  the  use  of  a digital  computer 
instead  of  the  analog  computer  highly  feasible. 

Background 

The  use  of  digital  flight  control  techniques  for  air- 
craft control  functions  has  been  shown  by  recent  research 
efforts  using  the  A7-D  and  F-8  aircraft (Ref  l;Ref  2).  There 
are  many  different  design  methods  for  digital  controllers. 
These  methods  fall  into  two  broad  categories:  (1)  digitiza- 
tion, which  is  the  discretization  of  a compensator  designed 
in  the  continuous  domain  (s-plane),  and  (2)  direct  digital 
design,  which  is  the  direct  design  of  the  compensator  in  the 
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discrete  domain  (z-  or  w-plane). 
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Digitization . In  a recent  report  by  Honeywell  (Ref  3), 
the  digitization  approach  using  the  Tustin  approximation  is 
advocated.  The  Tustin  transformation  is  a bilinear  transfor- 
mation technique  that  has  the  desirable  properties  of  cascade 
preservation,  stability  invariance,  and  steady-state  gain  in- 
variance. It  is  also  relatively  easy  to  apply  and  understand. 

Another  means  of  discretizing  an  s-plane  compensator  is 
the  z-transf ormation . This  method  has  the  same  properties  as 
previously  mentioned  for  the  Tustin  transformation  except  for 
cascade  preservation.  Of  the  numerous  digitization  techniques, 
the  Tustin  transformation  and  the  z-transf ormation  have  found 
the  widest  use. 


Direct  Digital  Design . As  in  the  digitization  process, 
there  are  several  methods  that  can  be  applied  to  accomplish 
a direct  digital  design.  These  methods  include,  but  are  not 
limited  to,  pole-placement  in  the  z- , w- , or  w' -plane  and  a 
discrete  design  using  optimal  techniques.  In  a report  by 
Powell,  the  discrete  method  was  applied  to  the  analysis  of  a 
digital  system  for  a Boeing  737  aircraft  with  favorable  re- 
sults (Ref  4). 

Scope 

The  scope  of  the  problem  addressed  by  this  study  is  the 
design  of  the  digital  control  laws  for  an  advanced  fighter 
aircraft  performing  the  terrain  following  task.  The  aircraft 
is  projected  to  be  somewhat  similar  to  the  F-15.  The  design 
will  be  carried  forth  using  the  discrete  optimal  method  es- 
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r-;/^  tablished  by  Katz  and  Powell  (Ref  5).  The  design  is  based  on 

reference  control  inputs  to  the  elevator  and  thrust  level  and 

reference  states  provided  by  a reference  command  generator. 

The  reference  command  generator  receives  its  input  from  a 

path  generator.  The  reference  command  follower  system  is 

I shown  in  Figure  1.  The  reference  states,  x^.,  and  reference 

commands,  u , are  assumed  to  be  available.  One  means  of  ob- 
. ’ r 

taining  these  values  is  shown  in  Chapter  V. 

The  problem  has  been  defined  for  the  terrain  following 

i 

mode.  This  mode  provides  maneuvering  in  the  vertical  plane 

only;  thus,  the  problem  is  limited  to  the  longitudinal  axis 

i 

of  the  aircraft.  The  path  tracking  problem  reduces  to  a reg- 
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ulator  problem  in  the  continuous  domain.  The  state  cost 
matrix,  g,  and  the  control  cost  mr.trix,  R,  can  then  be  de- 
termined for  the  quadratic  performance  index,  J.  With  the 
cost  values  established,  an  evaluation  of  performance  against 
C*  criteria  is  made,  where  C*  is  defined  as  a blend  of  pitch- 
rate  and  normal  acceleration  (Ref  6).  Using  the  cost  values 
from  the  continuous  system,  the  discrete  optimal  regulator 
problem  is  determined  with  sample  rate  as  a parameter.  The 
resulting  system  is  then  analyzed  on  a hybrid  computer.  The 
results  are  presented  using  s-plane,  z-plane,  and  time  response 
analyses  of  the  system. 

Assumptions 

(1)  The  aircraft  dynamics  can  be  satisfactorily  repre- 
sented by  the  linearization  of  the  aircraft  motion  about  a 
nominal  flight  condition. 

(2)  The  aircraft  stability  derivatives  are  valid  for 
the  flight  regime  investigated.  In  addition,  these  deriv- 
atives remain  essentially  the  same  along  the  reference  flight 
trajectory  as  for  the  nominal  flight  condition. 

(3)  The  optimal  gain  matrix  is  approximated  as  the  steady- 
state  solution  to  the  Riccati  equation. 

(4)  The  reference  states,  x^,  and  reference  commands, 
u^,  are  available. 

Thesis  Organization 

Chapter  II  presents  the  system  equations  in  the  form  of 
a linearized  perturbation  model.  Chapter  III  illustrates  the 
development  of  the  system  as  a continuous  regulator  problem 
used  in  determining  a satisfactory  performance  index.  In 
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Chapter  IV  the  discrete  deterministic  regulator  problem  is 
solved  at  various  sample  rates  using  the  continuous  perfor- 
mance index  developed  in  Chapter  III.  An  acceptable  sample 
rate  for  use  in  the  discrete  system  is,  thus,  obtained.  A 
development  of  the  necessary  equations  for  use  in  the  refer- 
ence command  generator  is  presented  in  Chapter  V along  with 
a hybrid  simulation  of  the  terrain  following  task.  Chapter 
VI  contains  the  conclusions  and  recommendations  brought  forth 
by  this  research  effort. 
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II . System  Equations 

The  reference  command  follower  system,  as  shown  pre- 
viously in  Figure  1,  can  be  reduced  to  a regulator  problem. 
This  is  done  by  setting  the  reference  command,  u^.,  and  ref- 
erence state,  x^,  equal  to  zero.  The  system  then  becomes  the 
system  as  shown  in  Figure  2.  Linearization  of  the  aircraft 
equations  of  motion  about  the  nominal  flight  condition  pro- 
vides a basis  for  the  calculation  of  the  feedback  gain  matrix, 
K.  The  state  space  representation  for  this  system  is  shown 
in  Figure  3,  wnere  the  open  loop  plant  can  be  expressed  as 


X = ^ + Bu 

a) 

2 = C X 

(2) 

For  the  purpose  of  this  study  it  is  required  that: 


Figure  3.  State  space  block  diagram  for  the  regulator  pro- 
blem. 

1.  All  states  are  available  for  measurement. 

2.  The  plant  and  control  laws  are  linear. 

3.  The  system  is  time  invariant. 

In  addition  the  following  techniques  are  used: 

1.  Only  the  infinite-time  regulator  is  considered. 

2.  The  performance  index  to  be  minimized  is 

J = + u^R  u;  dt  (.3) 

with  2 symmetric  non-negative  definite  and  R pos- 
itive definite. 

With  these  facts  stated  it  is  necessary  to  develop  the  re- 
quired system  equations. 

Linearized  Aircraft  Model 

The  Assumed  state  variable  vector  is 
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X = ^ 'u  w 


q 0 6e  AT 


U) 


wnere 

'u  (.ft/sec)  is  the  perturbation  in  airspeed, 
w (tt/sec)  is  the  perturbation  in  vertical  velocity, 
q (.rad/sec)  is  the  perturbation  in  pitch  rate. 

0 (rad)  is  the  perturbation  in  pitch  angle. 

6e  ^rad)  is  the  perturbation  in  elevator  deflection. 
at  (lbs)  is  the  perturbation  in  thrust, 
h (ft)  is  the  perturbation  in  distance  from  nominal 
flight  path. 

(ft/sec)  is  the  total  forward  nominal  velocity. 
Since  it  has  been  assumed  the  aerodynamic  coefficients  along 
the  reference  trajectory  remain  the  same  as  for  the  nominal 
flight  condition,  then 


6x  = 

15) 

6x  = 

A X 

+ B u 

(6) 

where 

5u  = 

(7) 

and 

u - 

5e 

c 

AT 

c 

• m 

^8) 

6e  is  the  commanded  elevator  deflection  in  radians  and  AT 
c c 

IS  the  commanded  thrust  change  in  pounds.  Non-zero  values 

for  X and  u are  considered  in  Chapter  V. 

— r ““r 
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Open  Loop  Dynamics  Matrix.  A,  and  Control  Distribution  Matrix . B 
In  order  to  solve  the  regulator  problem,  it  is  first  nec- 
essary to  obtain  the  system  equations  in  state  space  form  us- 
ing the  dimensional  stability  derivatives.  Tne  longitudinal 
perturbation  equations  of  motion  can  be  represented  by 

'u  = X^-u  + X^w  + X^q  - g 0 + X5g6e  + X^.j.AT 

w = Zu'"  V ^'^n  ^ 2q)  q + + Z^.j.AT 

q = M 'u  + Mw+Mq  + M-  6e 
^ u w q^  6e 

e = q . 

when  the  nominal  pitch  angle  is  zero  (Ref  7:  4.113  ).  The 
equation  for  the  elevator  servo  can  be  approximated  by  the 
transfer  function 


6e(s) 

6e^(sj 


1 


1 


which  can  be  written  as 


1 

6e  = 

^6e 


(10) 


(11) 


where  t.  is  the  servo  time  constant  and  6e  is  the  commanded 
6e  c 

change  in  elevator  deflection.  Similarly,  for  the  throttle 
servo 


AT(s) 

AT^is) 


1 

1 


a2) 


which  implies 
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AT  = ^ (AT^  - AT) 

^AT  ° 


wGere  is  the  thrust  servo  time  constant. 
The  relationship  for  h is 

h = U sin  0 - w cos  0 


(13) 


(14) 


where  0 = 0^  + 6 ,ie.,  nominal  pitch  attitude  plus  the  per- 
turoation  in  pitch.  For  tne  case  under  study,  0^  = 0.  thus 


h = U sin  0 - w cos  9 


(15) 


When  linearized,  equation  (15;  becomes 


h = V 0 - w 
n 


(16) 


Combining  equations  (9),  (11),  (13),  and  (16)  into  state 
space  form  yields 


X = A X + B u 


(1) 


and  the  control  distribution  matrix  is 


(18) 


Using  the  derivatives  provided  by  the  Avionics  System  Group, 
Aerospace  Systems  Division,  Wright-Patterson  Air  Force  Base, 
for  a flight  condition  of  .8  Mach  (528kts  or  882  ft/sec)  at 
sea  level,  the  open  loop  dynamics  matrix  and  control  distri- 
bution matrix  become  those  shown  in  Figure  4. 

Output  Matrix,  C 

The  given  outputs  of  the  system  are 


^ = C X 

which  for  the  regulator  problem  becomes 


(2) 


Z 

— 

C 6x 

(19) 

since 

6x 

s 

X 

(5) 

and 

iZ 

s 

Z 

(7) 

The  specific 

outputs 

of 

tne  system  are  defined  to  be 
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These  outputs  can  be  obtained  from  the  aircraft  navigation 
system  and  Central  Air  Data  Computer  (CAdC).  The  variables 
used  as- outputs  describe  the  aircraft  flight  path.  The  var- 
iable h is  the  deviation  from  the  nominal  altituae  as  pre- 
viously aefinea.  Tne  remaining  output  variables  are  defined 
as 


h ^ 

dh 

. dt 

dh 

dt 

dR 

dR 

(21) 


where  Y is  the  perturbation  in  flight  path  angle  or  slope 
of  the  flight  path,  and 


where  K is  the  flight  path  curvature,  and  R is  the  horizon- 
tal range.  It  is  preferable  to  work  in  range  intervals, 
rather  than  time,  for  tnen  the  usual  assumption  of  constant 
horizontal  velocity  is  not  required  to  relate  the  terrain 
locations  to  those  of  the  aircraft  (Ref  8;.  The  variable 
'u  IS  the  deviation  from  nominal  airspeed  as  previously  ae- 
f ined. 

With  these  definitions  in  mind,  it  is  now  possible  to 
determine  the  output  matrix,  C.  As  shown  in  equation  Cl6) 
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Y = 


w 

V 


+ e 


(24) 


Differentiating  with  respect  to  time  and  dividing  by  yields 


k = 


V* 

n 


w 

K 


(25) 


When  the  value  of  w from  equation  (9)  is  substituted  into 
equation  (25),  equation  (25)  becomes 


K 


■Zu'u  - - Z^q  - - Z^^AT 

''n 


(26) 


Equations  (24)  and  (26),  along  with  h and  'u  which  are  directly 
observable,  can  be  combined  to  form  the  output  matrix,  C. 


C = 


0 

0 


0 

V_ 


0 

0 


■Z  -Z  -Z 
u w 


V2 

n 


V2 

n 


V2 

n 


0 

1 


0 

0 


0 

0 


1 

0 


'^6e  "^AT 


V2 

n 


V2 

n 


0 


0 


(27) 


With  the  system  equations  so  defined,  it  is  now  possible  to 
proceed  with  the  solution  of  the  regulator  problem.  The  C 
matrix  for  the  aircraft  at  .8  Mach  and  at  sea  level  is  shown 
in  Figure  5. 
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III.  Continuous  Deterministic  Regulator  Problem 
The  regulator  problem  involves  a plant  whose  output,  or  ; 

I 

any  of  the  derivatives  of  the  output,  is  initially  nonzero. 

The  objective  is  to  provide  a control  to  take  the  plant  from  | 

this  nonzero  state  to  the  zero  state  and  to  provide  a satis-  1 

factory  transient  response.  This  is  accomplished  by  providing  j 

a control  composed  of  a linear  combination  of  the  states.  I 

Various  methods  have  been  presented  in  the  literature  for 

determining  the  necessary  feedback  gains  for  this  control.  For  ^ 

the  problem  under  study,  it  is  desired  to  minimize  a quadratic 
performance  index 

J = i X + u'^R  u)  dt  (28) 

Bryson  and  Hall  developed  a computer  program,  OPTSYS,  for 
the  solution  of  this  optimal  control  problem  using  eigenvector 
decomposition  (Ref  9) . This  method  was  reported  by  them  to 
be  faster  and  more  accurate  for  high  order  systems  than  the 
Ricatti  equation  method.  Thus,  the  OPTSYS  computer  program 

j is  suited  for  use  in  determining  the  feedback  gain  matrix,  K. 

I 

; The  control,  u,  for  the  regulator  becomes 

u = K X (29) 

where  u is  composed  of  the  change  in  elevator  command,  6e^, 
and  the  change  in  commanded  thrust  level,  AT^.  With  the 
values  for  the  feedback  gain  matrix  available  from  OPTSYS, 
given  a state  cost  matrix,  S,  and  a control  cost  matrix,  R, 
it  is  possible  to  obtain  root-lci  of  the  eigenvalues  as  a 
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function  of  the  elements  of  S and  R.  From  the  root-loci  and 
selective  analog  simulations,  it  becomes  possible  to  deter- 
mine a satisfactory  state  cost  matrix  and  control  cost  matrix 
for  use  in  the  discrete  optimal  control  problem. 

State  Cost  Matrix  Development 

For  the  terrain  folio”  g task,  it  is  desired  to  minimize 
the  energy  of  the  output,  2,  in  conjunction  with  the  control, 
u,  in  the  performance  index 

J = i /”  £ + u”^  R u)  dt  (3) 

Thus,  the  output  cost  matrix  of  equation  (3)  must  be  related 
to  the  state  cost  matrix  of  equation  (28)  prior  to  use  in  the 
OPTSYS  program.  Since 

2.  = C X (2) 

then  Z “ (£  S (£ 

= x'^  2 C X (30) 

and  the  state  cost  matrix  becomes 

S = 2 £ (31) 

Since  a diagonal  cost  matrix  is  often  satisfactory,  the  out- 
put cost  matrix  is  defined  to  be 
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(32) 


S = 


0 0 0 
0 q 0 0 


0 q 0 

3 


0 0 


Since  there  is  an  infinite  number  of  possible  choices 


for  the  q^  terms,  a means  of  judiciously  narrowing  these 
choices  must  be  used.  One  such  method  is  to  normalize  the 
maximum  RMS  value  of  the  output  variables  and  their  corres- 
ponding weighting  factors  to  unity  in  defining  the  values  for 
In  the  case  at  hand,  the  output  variables  are 


thus 


T 

Z Q Z 


• •• 

h S—  'u 

V 

n n 


= h^q.  . * f-!l-l=q.  'u 


(20) 


Y 2 3 

n 


(33) 


If  the  maximum  value  for  each  term  is  set  equal  to  one,  the 


q.  terms  can  be  defined  as 


m'  n 


(34) 


(35) 


(36) 


(37) 


where  the  subscript  m denotes  the  variables  maximum  value. 

It  can  be  observed  from  the  foregoing  that  eq^al  weight  has 
been  placed  on  the  maximum  RMS  value  of  each  output  variable. 
While  this  weighting  may  not  produce  a desirable  system,  it 
does  have  the  benefit  of  intuitive  appeal  and  does  yield  a 
starting  point  for  the  iterative  design  approach.  The  next 
step  is  to  select  what  is  felt  to  be  the  maximum  allow- 
able values  for  the  output  variables. 

The  maximum  deviation  from  commanded  flight  path,  h^, 
is  chosen  as  20  feet.  The  maximum  rate  of  change  in  the  de- 
viation  from  flight  path,  h^,  is  chosen  to  be  40  feet  per  se- 
cond.  For  the  acceleration  term,  h , a value  of  40  feet  per 
second  squared  is  chosen.  The  initial  value  for  the  maximiun 
deviation  in  airspeed,  u^,  is  100  feet  per  second.  These 
maximum  values  may  be  rather  large,  but  provide  a realistic 
starting  point  for  forming  the  root-loci  of  the  system. 

Using  the  values  presented  above,  the  output  cost  matrix 
becomes 


s = 


0.0025  0.0  0.0  0.0 


486.2  0.0  0.0 


0.0  3.78E08  0.0 


0.0  0.0  0.0001 


(38) 


The  state  cost  matrix  is  derived  from  the  relationship  shown 
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in  equation  (31).  For  the  values  of  £ in  equation  (38),  the 
state  cost  matrix  is  shown  in  Figure  6. 

Control  Cost  Mati  Development 

The  control  cost  matrix,  R,  is  initially  established  in 
the  same  manner  as  that  for  the  output  cost  matrix.  Namely, 


R = 


(39) 


where  r is  the  weight  on  the  control  energy  of  6e  and  r 

1 c 2 

is  the  weight  on  the  throttle  (thrust)  control,  AT^.  Thus, 


r 


1 


r 


2 


1 

(6e^  )2(Rad)2 
m 


(57.3) 

( 6e^  )2  (deg)2 


(40) 


1 


(AT  y 


(41) 


With  the  methods  established  for  determining  the  state  and 
control  cost  matrices,  it  is  now  possible  to  proceed  with  a 
root-locus  study  of  the  system. 

Root-Locus 

The  root-loci  for  the  system  are  determined  by  varying 
the  cost  functions  which  results  in  a variation  in  the  state 
feedback  gain  quantities  and,  thus  a change  in  the  location 
of  the  characteristic  roots.  In  the  following  figures  depict- 
ing the  root-loci,  a single  cost  parameter  is  varied  while 
the  other  parameters  are  held  constant.  The  associated  reg- 
ulator problem  for  each  variation  in  cost  is  then  solved  to 
yield  the  closed  loop  eigenvalues  (roots)  and  feedback  gains 
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using  the  OPTSYS  program.  The  root  locations  are  then  plotted 
as  a function  of  the  cost  parameter. 

Closed  Loop  Pole  Placement 

With  the  root-loci  as  design  tools,  the  closed  loop  roots 
can  be  placed  at  a desired  operating  point  by  a judicious  choice 
of  the  state  and  control  cost  parameters.  Using  the  Cornell 
Aeronautical  Laboratory  "thumbprint"  shown  in  Figure  7 as  a 
guide,  a damping  ratio  of  approximately  .7  and  a natural  fre- 
quency of  five  radians  per  second  was  chosen  as  the  desired 
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Figure  8.  Open  loop  eigenvalues  and  eigenvectors. 


operating  point  for  the  short  period  roots.  The  open  loop 
eigensystem  is  depicted  in  Figure  8.  The  root-loci  are  shown 
in  Figures  9 through  14. 


As  shown  in  Figure  9,  if  the  cost  parameter  is  set 
equal  to  6.72E+10,  h^  equal  to  three  feet  per  second  squared, 
the  desired  root  location  has  been  approximated.  This  root 
location  is  -3.8  ± j 3.44  with  a damping  ratio  of  .74  and  a 
natural  frequency  of  5.1  radians  per  second.  As  indicated  in 
Figure  10,  the  control  cost  parameter  r^  could  also  be  used 
either  independently  or  inconjunction  with  q to  achieve  the 

3 

desired  pole  placement.  The  problem  encountered  in  using 
the  control  cost  parameter  is  that  to  achieve  the  desired  pole 
placement  would  require  additional  control  energy.  It  is  felt 
at  this  point  that  the  one  and  one-half  degree  elevator  deflec- 
tion, as  it  relates  to  energy,  is  sufficient. 

Figure  11  shows  the  root-locus  as  a function  of  state 

cost  parameter  q , the  weight  on  the  h term  of  the  output. 

1 

A decrease  in  the  q^  term  can  be  thought  of  in  direct  cor- 
relation to  how  fast  the  regulator  drives  an  initial  condi- 
tion on  the  state  h toward  zero.  Initially,  the  phugoid  root 
will  be  placed  at  -.273  ± j .268  which  has  a damping  ratio 
of  .71  and  a natural  frequency  of  .38  radians  per  second. 

This  corresponds  to  q^  equal  to  2.5E-03  which  is  based  on 
hjjj  equal  to  twenty  feet  as  derived  previously  in  this  chap- 


ter. 


The  remaining  roots  of  interest  are  the  roots  that  cor- 


respond to  the  airspeed  and  thrust  response.  As  shown  in 
Figures  12  and  13,  these  root  locations  can  be  varied  using 
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Figure  9.  Root-locus  as  a function  of  state  cost  qj  with 
Qi  = .0025,  q2  = 486.2,  q4  = .0001,  Tj  = 1459, 
rj  = 6.25E-08. 


Figure  10.  Root-locus  as  a function  of  control  cost  rj  with 

qi  = .0025,  qz  = 486.2,  qj  = 3.78E+08,  q^  = .0001, 
T2  = 6.25E-08. 
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Figure  11.  Root- locus  as  a function  of  state  cost  Qj  with 

Qz  = 486.2,  Qs  = 3.78E+08,  = .0001,  Tj  = 1459, 

r2  = 6.25E-08. 


o Open  Loop 
X Closed  Loop 


m 


Figure  12.  Root-locus  as  a function  of  state  cost  q4  with 

qi  = .0025,  qa  = 486.2,  qs  = 3.78E+08,  rj  = 1459, 
ra  = 6.25E-08. 
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Figure  14.  Root-locus  as  a function  of  state  cost  Qj  with 

Qj  = .0025,  Qj  = 3.78E+08,  q,,  = .0001,  Tj  = 1459, 
Pj  = 6.25E-08. 
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either  the  parameter  or  r^.  Since  for  this  study  it  is 
desired  to  hold  throttle  cycling  to  a minimum,  the  closed 
loop  location  of  the  roots  that  correlate  with  the  airspeed 
and  thrust  response  will  be  established  by  the  original  value 
of  q equal  to  .0001.  The  parameter  r will  be  set  equal  to 

4 2 

6.25E-08  which  corresponds  to  4000  pounds  for  AT^.  The  4000 
pound  figure  is  taken  as  ten  per  cent  of  the  available  mil- 
itary thrust  at  sea-level.  Figure  15  shows  the  output  cost 
matrix,  the  corresponding  state  cost  matrix,  S,  the  con- 
trol cost  matrix,  R,  and  the  feedback  gain  matrix,  K,  for 
the  above  root  locations.  Figure  16  lists  the  eigensystem 
of  the  optimal  closed  loop  system. 

Analog  Computer  Simulation 

The  linearized  longitudinal  equations  of  motion  devel- 
oped in  Chapter  II  were  programmed  on  an  analog  computer. 

The  basic  aircraft  response  to  a one  degree  elevator  com- 
mand of  one  second  duration  is  shown  in  Figure  17.  The  re- 
sponse graphically  displays  the  underdamped  nature  of  the 
short-period  response  (damping  ratio  equal  to  .22).  The 
control  law  was  then  implemented  using  the  feedback  gain 
matrix,  K,  shown  in  Figure  15. 

The  closed  loop  system  was  evaluated  for  an  Initial 
condition  on  altitude  deviation,  h,  airspeed  deviation, 

'u,  and  pitch  deviation,  6.  These  initial  conditions  were 
applied  individually  and  the  response  of  the  states  ob- 
served to  determine  the  "goodness"  of  the  control  law.  As 
shown  in  Figures  18,  19,  and  20,  the  state  responses  are 
well  behaved  for  each  of  the  initial  conditions.  As  a side 
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Figure  15.  Output  cost  matrix,  state  cost  matrix,  control  cost  matrix,  and  feedback  gain 
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Figure  16.  Closed  loop  eigensystem. 
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6e  (degrees)  9 (degrees)  q (degrees/sec)  w (ft/sec)  'u  (ft/sec) 


System  closed  loop  response  to  an  initial  condition 
of  20  feet  on  the  state  variable  h. 
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Figure  20.  System  closed  loop  response  to  an  initial  condition 
of  two  degrees  on  the  state  variable  6. 
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note  of  controller  "goodness",  the  controller  was  modified 
to  generate  a C*  response.  To  do  this  the  following  devel- 
opment was  necessary. 

C*  Time  History 

The  C*  criteria  was  first  proposed  by  the  Boeing  Air- 
craft Company  (Ref  6)  and  was  also  explored  as  a handling 
qualities  criteria  in  the  Survivable  Flight  Control  System 
study  conducted  in  1971  (Ref  10).  C*  is  defined  as 

C*  = kia^j  + k2e  (42) 


where  ki  and  k?  are  constants  and  a is  normal  acceleration, 

n 

Normal  acceleration  can  be  approximated  by 


Y 

a = n 
n 

32.2 


(43) 


where  y is  the  time  rate  of  change  of  the  flight  path  angle. 


Also 


thus 


32.2 


(44) 

(45) 


From  equation  (26) 


h = -Z'u-Zw-Zq-  Zr^6e  - Z._AT  (46) 

u w q oe  AT 


therefore 


C*  = 


ki(h) 


32.2 


+ k2q 


(47) 


The  steady  state  relationship  between  q and  a^  is 
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(32.2) 


(48) 


nss 


‘ss 


If  equal  importance  is  placed  on  the  steady  state  normal  ac- 
celeration and  pitch  rate  at  a mid-range  velocity  such  that 


k la 


nss 


kjq 


ss 


(49) 


then  the  mid-range  velocity  is  defined  as  the  crossover  vel- 
ocity, Solving  for  ki  yields 

k2(32.2) 

ki  = (50) 


and  C*  becomes 


C* 


(h)  + kjQ 


(51) 


Defining  a scaled  C*  response  as 

C* 


^5 


(V  ) 
' co' 


(52) 


then 

c.  = -Z^'U  - Z^w  f - Z^)q  - Zj^Se  - Z^^OT  (53) 

It  is  now  possible  to  obtain  an  equivalent  C*  time  history 
for  the  system  based  on  a selection  of  the  crossover  velocity. 

The  control  law  for  the  system  must  first  be  modified 
so  as  not  topenalize  for  a deviation  in  altitude  or  pitch. 

This  is  accomplished  by  setting  the  value  of  the  feedback 
gain  on  these  two  states  to  zero.  In  addition,  the  C*  cri- 
teria is  based  on  a short  period  approximation;  thus,  the 
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Figure  21. 


Normalized  C*  response  with  q2  = 486.2,  qj  = 6.72 
E+10,  qi,  = .0001,  ri  = 1459,  rz  = 6.23E-06  for 
two  values  of  cost  parameter  q . 


deviation  in  airspeed  is  defined  to  be  zero.  With  these  al- 
terations made  in  the  analog  program,  the  C*  time  history  is 
shown  in  Figure  21  with  C*  normalized  to  unity  for  the  steady 
state.  The  system  was  connected  in  the  command  follower  con- 
figuration and  a step  elevator  command  applied.  The  cross- 
over velocity  was  set  equal  to  400  feet  per  second. 

As  shown  in  Figure  21,  the  rise  time  for  the  C*  response 
is  closely  associated  with  the  cost  parameter  qj.  While  the 
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cost  matrices  were  chosen  to  provide  a desirable  state  re- 
sponse for  the  regulator  problem,  it  is  interesting  to  note 
the  "goodness"  of  the  selection  when  evaluated  against  the 
C*  criteria. 

Summary 

The  continuous  deterministic  regulator  problem  was 
solved  to  provide  the  state  cost  matrix  and  control  cost 
matrix  for  the  desired  closed  loop  system.  The  selection 
of  these  matrices  was  based  on  the  Cornell  Aeronautical  Lab- 
oratory "thumbprint"  which  specified  the  short  period  natural 
frequency  and  damping  ratio.  Once  the  selection  of  the  matrices 
was  completed,  an  analog  computer  simulation  was  accomplished 
to  verify  the  "goodness"  of  the  selection.  The  control  law 
was  then  modified  to  yield  a €♦  response  which  was  compared 
with  the  C*  criteria  developed  in  the  Survivable  Flight  Con- 
trol System  study  (Ref  10).  With  the  cost  matrices  deter- 
mined, it  is  now  possible  to  proceed  with  the  discrete  optimal 
regulator  problem. 


IV.  Discrete  Deterministic  Regulator  Problem 

There  are  many  different  methods  that  can  be  used  in  de- 
signing a discrete  compensator  for  a continuous  system  as  ex- 
pained  in  Chapter  I.  They  are  normally  broken  down  into  two 
categories:  (1)  digitization  of  a compensator  designed  in  the 
continuous  domain  (s-plane);  or,  (2)  the  direct  design  of  the 
compensator  in  the  discrete  domain  (z-  or  w-plane).  The  di- 
rect digital  design  state  space  approach  has  the  advantages 
of  not  introducing  errors  due  to  approximating  the  zero  order 
hold  as  in  the  digitization  procedure  and  the  state  space  ap- 
proach easily  handles  multi-input  multi-output  systems. 

In  this  chapter  the  discrete  deterministic  regulator 
problem  will  be  solved  as  a function  of  sample  rate.  The 
minimum  acceptable  sample  rate  will  be  based  on  the  effect 
the  sample  rate  has  on  the  closed  loop  roots.  Katz  and  Powell 
developed  a computer  program,  DISCUS,  for  the  solution  of  the 
discrete  regulator  problem  (Ref  5).  The  calculation  of  the 
regulator  is  based  on  eigenvector  decomposition  of  the  re- 
lated state-costate  Hamiltonian  matrix  similar  to  the  earlier 
work  by  Bryson  and  Hall  for  the  continuous  problem  (Ref  9). 

Included  in  the  DISCUS  program  are  subroutines  for  discretiz- 
ing a continuous  performance  index  and  an  associated  contin-  i 

uous  plant.  Thus,  the  state  cost  and  control  cost  matrices 
developed  in  Chapter  III  can  be  used  in  the  discrete  problem. 

While  it  is  possible  to  define  the  cost  matrices  directly  in  j 

the  discrete  domain  by  the  use  of  DISCUS,  the  insight  afforded  i 

by  the  continuous  domain  would  be  lost. 
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Discretization  of  the  Continuous  Plant 

The  discrete  formulation  of  a linear,  continuous,  dif- 
ferential system  can  be  thought  of  as  the  solution  in  the 
time  domain  from  one  point  in  time  to  another.  The  state 
model  of  a continuous  system  is  given  by  (repeated  from  Chap- 
ter I) 


X = A X + B u 


C X 


(1) 


(2) 


The  solution  to  equation  (1)  is 

X = exp  A(t  - t^)  x^(t^)  + exp  A(t  - t)  B u(T)dT 

— ~ O O 

(54) 

Assuming  that  the  inputs  are  held  constant  over  intervals  of 
T seconds,  the  input  u(t)  is  of  the  form 


u(t)  = u(kT),  kT  < t < (k+l)T 


(55) 


The  initial  condition  at  the  begining  of  the  interval  (kT, 
i*(k+l)T}  is 


xC^q)  = x(kT) 


To  determine  x(t)  at  the  end  of  the  interval 


t = kT 
o 


(56) 


t = (k+l)T 


is  substituted  into  equation  (54).  Thus 


(57) 

(58) 


(k+l)T 

X {(k+l)T}  = exp  ATx(kT)  +/  exp  A{(k+1)T  - t}B  u(kT)  dT 

■"  kT  (59) 
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Since  the  input  u(kT)  is  constant  for  the  integration  interval 
and  the  integration  is  valid  for  all  k,  by  making  the  change 
of  variable  t = (k+l)T  - T equation  (59)  becomes 

X {(k+l)T}  = exp  AT  x(kT)  + f^exp  At  B dt  ’ u(kT) 

« (60) 

Provided  the  sampling  interval  T is  constant  and  the  A and 
B matrices  are  time  invariant,  then 


exp  At 

= ±(  t) 

0 ^ T ^ T 

(61) 

and 

f^exp  At  B dt 

0 “ “ 

= G(t) 

(62) 

When  T = T,  ^ and  G are  constant  matrices  whose  elements  are 
functions  of  the  sampling  interval  T.  In  conclusion,  the 
discrete  state  model  for  a continuous  system  with  a piece- 
wise  constant  input  is 


X {(k+l)T} 

= ±(T)  x(kT)  + G(T)  u(kT) 

(63) 

Z(kT) 

= C x(kT) 

(64) 

Discretization  of  a Continuous  Performance  Index 
A continuous  quadratic  performance  index 


J 


X + u'^R  u)  dt 


(65) 


can  be  transformed  into  a discrete  version 


(66) 
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where  N = t/T  by  the  following  manipulations.  Equation  (65) 
can  be  written  as 


N-1  (i+l)T  T,  T, 

J = I / (X  S X + u^R  u)  dr  (67) 

i=o  iT 


since  x(t) 


_£(T-iT)  x^  + G(T-iT)  u^ 


(68) 


where  iT  t £ (i+l)T 

then 

(i+l)T  „ T 

/ S X + u R u)  dt  = 

iT 


All 

Al  2 

1 

•H 

XI 

1 

L ^ y 

Az  1 

CM 

<1 

— 1 
•H 

(69) 


where 


Ai  1 
Az  2 
Al  2 
Az  1 


T 

/ i'^(t)  S 1(t)  dr  (70) 


T 

S G(t)  + R}dT  (71) 
T 

/ 4>'^(t)  S G(t)  dT  (72) 

^o  — — — 

aTz  (73) 


The  DISCUS  program  has  subroutines  based  on  the  above  devel- 
opments for  discretizing  a continuous  plant  and  a continuous 
performance  index.  The  closed  loop  roots  and  discrete  feed- 
back gain  matrix  for  the  regulator  problem  can  now  be  obtained. 
Discrete  System  Closed  Loop  Roots 

Using  the  DISCUS  program,  the  discrete  regulator  prob- 
lem was  solved  at  various  sample  rates.  The  continuous  state 
and  control  cost  matrices  of  Chapter  III,  Figure  15,  were 
used  in  the  performance  index  along  with  the  sample  rates  of 
interest.  The  resulting  closed  loop  root  movement  for  the 
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Figure  22.  Discrete  closed  loop  root  movement  as  a function 
of  sample  rate. 


discrete  roots  (z-plane)  is  shown  in  Figure  22.  The  equiva- 
lent continuous  root  movement  (s-plane)  is  shown  in  Figure  23. 


Z-Plane  Roots . As  illustrated  in  Figure  22,  the  discrete 
closed  loop  root  movement  for  the  short  period  and  phugoid 


modes  forms  loci  with  a damping  ratio  of  approximately  .7. 
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Figure  23,  Equivalent  s-plane  discrete  closed  loop  root  lo- 
cation as  a function  of  sample  rate. 


This  is  true  of  the  short  period  root  movement  down  to  a fre- 
quency of  five  hertz.  This  should  be  the  case,  since  the  dom- 
inant roots  have  a damped  natural  frequency  of  3.44  radians 
per  second  (.55  Hz).  As  long  as  the  system  is  being  sampled 
at  a rate  at  least  ten  times  the  dominant  mode,  the  effects 
of  the  sampling  should  be  minimal.  When  the  sample  rate  is 
reduced  below  five  hertz,  the  locus  diverges  more  and  more 
from  the  .7  damping  ratio  indicating  the  effects  of  sampling 
at  less  than  ten  times  the  dominant  mode.  For  the  system  un- 
der study,  it  has  been  assumed  that  the  aircraft  is  a rigid 
body;  thus,  the  body  bending  modes  have  been  ignored.  The 
body  bending  modes  may  be  such  that  the  five  hertz  sample  rate 
would  not  excite  them.  If  this  is  not  the  case  the  rate  of 
sampling  would  have  to  be  adjusted.  This  particular  facet  of 
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of  the  problem  will  not  be  explored  further  by  this  study. 


Equivalent  S-Plane  Roots . The  equivalent  s-plane  roots 
are  determined  by  the  relationship 

sT 

z = e®^  (74) 

thus  s = (In  z)/T  (75) 

where  T is  the  sampling  period  and  s is  the  Laplace  operator 
The  periodic  nature  of  equation  (75)  must  be  realized  to 
fully  describe  the  equivalent  s-plane  roots.  Since  for  a 
complex  variable 

In  z = In  |z|  + j arg  z ± 2nirj 

(76) 

for  n = 0,  1,  2,  

the  complex  natural  logarithm  is  infinitely  many  valued. 
Figure  23  indicates  only  the  portion  of  the  primary  strip 
that  contains  the  principle  value  of  the  logarithm.  It  is 
pointed  out  that  to  use  an  equivalent  s-plane  mapping,  the 
effects  of  the  folded  roots,  the  roots  indicated  by  the 
±2nTTj  term  in  the  complex  logarithm,  must  also  be  considered 
At  a sample  rate  of  five  hertz,  the  primary  strip  extends 
from  -15.7  to  1-15.7  radians  per  second.  Thus,  the  nearest 
folded  root  is  located  at  -3.8  ± j 28.  This  is  far  enough 
removed  from  the  dominant  roots  to  have  minimal  effect  on 
the  system.  This  is  not  true  of  a two  hertz  sample  rate 
since  the  primary  strip  has  moved  into  ±6.28  radians  per 
second.  At  two  hertz  the  nearest  folded  roots  are  located 
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at  -3.8  ± j 9.1  and  may  degrade  the  response  of  the  system. 

It  is  also  shown  in  Figure  23  that  at  five  hertz  the  system 
closed  loop  root  location  closely  approximates  the  continuous 
root  location. 

Hybrid  Regulator  Simulation 

The  linear  perturbation  model  developed  in  Chapter  II 
was  programmed  on  an  analog  computer  and  the  control  laws 
were  implemented  on  a Xerox  Sigma  7,  112K,  32-bit  word  com- 
puter. The  systems  were  interfaced  using  +100  volt,  15  bit, 

analog  to  digital  converters  and  ±100  volt,  15  bit  digital 

« 

to  analog  converters. 

The  closed  loop  system  in  the  form  of  the  discrete  reg- 
ulator was  evaluated  for  an  initial  condition  on  altitude  de- 
viation, h,  airspeed  deviation,  'u,  and  pitch  deviation,  0. 

The  initial  conditions  were  applied  individually,  as  in  the 
evaluation  of  the  continuous  regulator  of  Chapter  III,  using 
sample  rates  from  100  hertz  to  2 hertz.  The  states  were  ob- 
served to  determine  the  effects  of  the  sampling  in  the  time 
response.  As  expected,  at  the  higher  sample  rates  no  discern- 
able  difference  from  the  continuous  system  could  be  observed. 
As  the  sample  rate  was  lowered  below  20  hertz,  the  discrete 
nature  of  the  control  became  apparent.  A slight  ripple  was 
observable  in  some  of  the  states.  While  this  ripple,  or  noise, 
was  not  predictable  on  the  basis  of  the  z-  and  s-plane  root 
analysis,  it  is  pointed  out  that  these  analysis  tools  do  not 
take  into  account  the  finite  word  length  constraint  of  the 
digital  computer  and  the  associated  interfacing.  While  this 


computer  had  a word  length  of  32  bits,  the  relatively  small 


word  length  of  the  analog  to  digital  and  digital  to  analog 
converters  introduced  a round  off  error  that  was  evident  at 


the  lower  sample  rates.  At  the  higher  rates  of  sampling,  the 
noise  was  effectively  filtered  by  the  elevator  and  throttle 
actuators . 

Figure  24  shows  the  system  state  response  at  five  hertz 
to  an  initial  condition  of  20  feet  on  the  state  variable  h. 

The  noise  introduced  by  the  discrete  system  is  most  evident 
in  the  state  variable  AT.  The  noise  appears  to  have  a maxi- 
mum strength  of  plus  or  minus  one  pound  - which  is  small  in- 
deed. The  regulator  response  to  an  initial  condition  on  the 
state  variable  'u  is  shown  in  Figure  25.  The  noise  is  no 
longer  apparent  in  the  AT  state  due  to  the  scaling  of  the  var- 
iable for  recording.  The  noise  is  noticeable  in  e,  q,  and 
fie.  The  noise  has  a maximum  strength  of  approximately  ±.01 
degree  in  the  elevator  deflection,  ±.05  degree  in  pitch,  and 
±.01  degree  per  second  in  pitch  rate.  The  commanded  elevator 
deflection  is  included  in  Figure  25  to  illustrate  the  discrete 
nature  of  the  noise  introduced  by  the  digital  process.  Figure 
26  depicts  the  discrete  regulator  response  at  a sample  rate  of 
five  hertz  to  an  initial  condition  of  two  degrees  on  the  state 
variable  9.  The  noise  is  suppressed  in  this  figure  due  to  the 
lower  amplifications  required  in  scaling  the  variables  for  re- 
cording. 

The  state  variable  closed  loop  responses  at  a sample  rate 
of  two  hertz  for  the  same  initial  conditions  are  shown  in  Fig- 
ures 27,  28,  and  29.  The  magnitude  of  the  noise  has  increased 
since  the  error  introduced  in  the  control  law  has  longer  to 
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act  in  the  system  before  the  system  is  again  sampled. 
Summary 


Direct  digital  design  has  the  advantage  of  not  intro- 
ducing errors  into  a discrete  system  since  the  zero  order 
hold  is  accounted  for  precisely.  The  state  space  approach 
lends  itself  nicely  to  multiple  input,  multiple  output  sys- 
tems and,  through  the  use  of  optimal  control  techniques,  the 
design  engineer  can  readily  come  to  grips  with  very  demanding 
control  problems.  It  was  shown  in  this  chapter  how  an  opti- 
mal regulator  problem  was  discretized.  The  effects  of  sample 
rate  on  the  system  were  then  examined.  The  finite  word  length 
of  a digital  system  was  shown  to  produce  the  equivalent  of  a 
noise  being  introduced  into  the  system.  The  sampling  rate  of 
five  hertz  was  shown  to  be  high  enough  to  form  the  desired 

( 

controls.  In  the  following  chapter,  the  task  of  terrain  fol-  i 

lowing  will  be  explored  using  the  regulator  that  has  thus  far  | 

been  developed. 


f 

I 

i 
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Figure  24.  System  closed  loop  response  to  an  initial  condition 
of  20  feet  on  the  state  variable  h using  a five 
hertz  sample  rate. 
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(degrees)  (degrees)  (pounds)  (degrees)  (deg/sec) 


System  closed  loop  response  to  an  initial  condition 
of  -20  feet  per  second  on  the  state  variable  'u 
using  a five  hertz  sample  rate. 


(feet)  I (pounds)  (degrees)  (degrees)  (deg/sec) 


System  closed  loop  response  to  an  initial  condition 
of  two  degrees  on  the  state  variable  9 using  a sample 
rate  of  five  hertz. 


(feet)  (pounds)  (degrees)  (degrees) 


System  closed  loop  response  to  an  initial  condition 
of  20  feet  on  the  state  variable  h using  a two  hertz 
sample  rate. 


Figure  28.  System  closed  loop  response  to  an  initial  condition 
of  -20  feet  per  second  on  the  state  variable  'u  us- 
ing a two  hertz  sample  rate. 


(degrees)  (degrees) 


Figure  29.  System  closed  loop  response  to  an  initial  condition 
of  two  degrees  on  the  state  variable  6 using  a sam- 
ple rate  of  two  hertz. 
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V.  The  Terrain  Following  Task 

In  addition  to  the  regulator  developed  in  the  previous 
chapters,  it  is  necessary  to  derive  the  reference  command 
generator  and  to  obtain  the  desired  path  information  in  or- 
der to  implement  the  system  for  the  terrain  following  task. 

For  the  purpose  of  this  study,  the  desired  path  over  the 
terrain  is  assumed  to  be  generated  by  an  optimal  cubic  spline 
algorithm  (Ref  11).  The  algorithm  stores  the  desired  alti- 
tude and  the  first  derivative  of  altitude  with  respect  to 
horizontal  range  indexed  at  prescribed  range  intervals.  The 
higher  derivatives  of  the  path  could  also  be  stored  but  at 
an  additional  cost  in  required  memory.  This  chapter  will  il-  I 

lustrate  one  method  of  reconstructing  the  desired  path  and 
path  derivatives.  Using  the  reconstructed  path  information, 
the  reference  commands  and  reference  states  are  generated. 

1 

A change  of  notation  is  required  since  the  total  states  j 

I 

and  commands  are  the  quantities  of  interest.  From  this  point  | 

on  in  the  thesis,  the  notation  will  refer  to  the  total  var-  t 

iable  instead  of  the  perturbation  about  the  nominal  condition 
as  used  previously.  In  the  later  portion  of  the  chapter,  the 
total  system  is  evaluated  through  the  use  of  a hybrid  com- 
puter. A development  of  the  non-linear  aircraft  longitudinal 
equations  of  motion  is  included  in  Appendix  A.  These  equations 
were  programmed  on  the  analog  portion  of  the  hybrid  system  to 
serve  as  the  truth  model  for  the  system  under  study.  Included 
in  Appendix  B is  a computer  program  REFCOM,  developed  to  re- 
construct the  path  and  furnish  the  commands  and  reference 
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states.  The  feedback  gain  quantities  for  five  of  the  sample 
rates  employed  are  listed  in  Appendix  C. 

Path  Reconstruction 


I 

i 


t 


It  is  shown  in  the  literature  (Ref  11;  152-157)  how  a 
cubic  spline  curve  may  be  reconstructed  and  that  derivation 
will  not  be  repeated  here.  The  required  difference  equations 
to  perform  this  reconstruction  from  a known  altitude,  slope, 
and  range  are 


h^(a) 

h*(o) 

h-'(a) 


+ 


+ 


+ 


h(k)  + a^(3  - 2a){h(k+l)  - h(k)} 

{(a  - 20^  + 03)h'(k)  + (a^  - a2)h'(k+l)}AR  (77) 

-^(1  - a){h(k+l)  - h(k)}  + {(1  - 4a  + 3a2)h’(k) 
AR 

(30^  - 2o)h’(k+l)}  (78) 

■ 6(1.=  2ai  {ii(k+l)  - h(k)}  + {(3o  - 2)h'(k) 

(3a  - l)h’(k+l)}  (79) 


h^’ '(a) 


-12 
(AR)  = 


{h(k+l)  - h(k)}  + 


(AR) 


{h'(k) 


+ h'(k+l)} 


(80) 


where  a is  the  non-dimensional  range  variable  defined  by 


R - R(k) 
AR 


(81) 


The  variable  R is  the  current  range  and  AR  is  the  dimensional 
range  increment  expressed  by 
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r 


AR 


R(k+1)  - R(k) 


(82) 


The  primes  in  equations  (77)  through  (80)  indicate  the  first 
through  third  derivatives  of  altititude,  h,  with  respect  to 
horizontal  range.  The  parenthetical  terms  (k)  and  (k+1)  de- 
note the  current  interval  of  operation  as  defined  by  the 
range  variable  R.  Thus,  the  path  is  completely  reconstructed 
for  any  range  given  the  discrete  range  points  R(k)  and  R(k+1), 
along  with  the  slope  at  these  points,  h'(k)  and  h'(k+l),  and 
the  path  altitude  at  these  points,  h(k)  and  h(k+l).  It  is 
now  necessary  to  define  the  reference  command  generator. 
Reference  Command  Generator 

The  reference  command  generator  computes  the  desired 
input  commands  and  reference  states  for  use  in  the  command 
follower  system  of  Chapter  I.  A simple  form  of  constant  energy 
control  is  applied  to  determine  the  desired  forward  velocity, 
V^.  Since  the  energy  of  the  system  is  held  constant  over  the 
path  legs,  the  resulting  control  system  does  not  require  a 
constant  velocity.  Throttle  cycling  is  minimized  and  engine 
life  increased.  The  desired  energy  level  for  each  leg  could 
be  either  calculated  using  apriori  knowledge  of  the  terrain 
or  updated  for  each  path  leg  based  on  the  terrain  of  the  pre- 
vious leg  and  the  required  time  over  the  path  end  point. 

The  length  of  each  path  leg  would  be  determined  by  the  nature 
of  the  mission  and  the  required  tolerance  in  the  time  over 
target.  If,  due  to  navigational  restraints,  a constant  air- 
speed controller  were  desired,  the  reference  generator  could 
be  greatly  simplified.  The  following  is  one  means  of  deriving 
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a constant  energy  reference  command  generator. 


Desired  Velocity . The  energy,  E,  of  the  center  of  mass 
required  for  a path  leg  can  be  defined  by  the  sum  of  the  kin- 
etic and  potential  energies.  The  following  represents  this 
relationship 

E = J mV*  + mgh  (83) 

If  mass,  m,  is  assumed  constant,  a scaled  version  of  total 
energy  can  be  obtained.  Dividing  both  sides  of  the  energy 
equation  by  a scaling  factor  equal  to  mass,  the  energy  re- 
lationship becomes 

i ''o  * 

where  V^  is  the  average  forward  velocity  required  to  achieve 
the  time  over  the  path  leg  end  point  and  h^  is  mean  path  al- 
titude. The  average  forward  velocity  is  equivalent  to  the 
nominal  velocity,  V^,  used  in  the  derivation  of  the  feed- 
back gain  quantities  in  the  previous  chapters.  The  rela- 
tionship for  the  desired  velocity  becomes 

Vj  = V2(Eo  - el'd’ 

where  h^  is  the  desired  path  altitude  provided  by  the  path 
generator.  The  desired  velocity  and  its  associated  deriva- 
tives are  required  in  the  development  of  the  reference 
states  and  controls.  The  d subscript  indicating  a desired 
quantity  will  be  dropped  from  this  point  on  in  the  derivation 

except  where  needed  for  clarity.  Since  E is  constant,  the 
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first  derivative  of  desired  velocity  is 

gh 

V 

V 


(86) 


The  path  generator  furnishes  the  spacial  derivatives  of  al- 
titude and  these  derivatives  are  related  to  the  time  deriv- 
ative by 

dh  dR 

h = • 

dR  dt 

= h'  V cos  Y (87) 

dh 

where  h ’ = 

dR  (88) 

dR 

and  = V cos  Y 

dt 

The  term  Y is  the  desired  flight  path  angle.  Substituting 
the  above  value  for  h into  equation  (86),  V becomes 


V 

= - g h'  cos 

Y 

(90) 

Since 

h’ 

= tan  Y 

(91) 

an  equivalent 

expression 

for  the  derivative  of 

desired  vel- 

ocity  is 

V 

= - g sin  Y 

(92) 

Taking  the  second  derivative  yields 

V 

= - g Y cos 

Y 

(93) 

Since 

Y 

= arc  tan  h' 

(94) 

59 


r 1 


1 dh' 

then  y (95) 


1 + (h' )2  dt 

Using  the  further  relationships 

sin^Y 


(h-  )2 

= 

tan^Y  = 

(96) 

cos^Y 

dh' 

and 

dt 

h' ' V cos  Y 

(97) 

Y becomes 

y 

= 

h ' ' V cos  *Y 

(98) 

Defining 

K 

= 

h ' ' cos  ®Y 

(99) 

then 

Y 

= 

K V 

(100) 

And  the  second  derivative 

of  the  flight  path  angle 

is 

Y 

= 

K V + K V 

(101) 

where 

K 

= 

V(h' ' ' cos^Y  - SK^h' ) 

(102) 

Substituting 

the  value 

of 

Y into  equation  (93),  V 

becomes 

V 

= 

- g K V cos  Y 

(103) 

The  desired  velocity  and  its  derivatives  have,  thus,  been 
defined  in  terms  of  the  available  quantities  from  the  path 
generator.  It  is  next  necessary  to  derive  the  reference 
angle  of  attack,  a^.,  and  elevator  reference  command,  60^.. 

Once  this  is  accomplished,  the  remaining  elements  of  the 
state  reference  vector  can  be  derived. 

Reference  Angle  of  Attack  and  Elevator  Command . In  order 
to  solve  for  the  reference  angle  of  attack  and  elevator  coro- 
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mand,  it  is  necessary  to  formulate  two  equations  with  these 
two  quantities  as  unknowns  for  simultaneous  solution.  This 
can  be  done  by  examining  the  aircraft  equations  of  motion  in 
conjunction  with  the  equations  for  aircraft  lift  and  moment. 

The  sum  of  the  forces  in  the  direction  perpendicular  to 
the  free  stream  velocity  can  be  expressed  by 


m a = L - mg  cos  y.  + T sin(e  + a ) (104) 

n Cl  6 z* 

where  a is  normal  acceleration,  L is  lift,  T is  the  thrust 
estimate  (this  quantity  will  be  approximated  by  the  thrust 
required  from  the  previous  instant),  and  e is  the  thrust  line 
offset.  Using  a small  angle  approximation  for  the  sine  term 
simplifies  the  equation  and  is  a valid  assumption  provided 
the  angle  of  attack  and  thrust  line  offset  are  small.  For 
a "hard"  ride  during  the  terrain  following  task,  the  optimal 
spline  path  is  restricted  to  -1  "g"  and  plus  2 "g"  from 
nominal  (Ref  11).  With  this  restriction,  the  angle  of  attack 
is  limited  to  a small  angle  validating  the  assumption. 

The  lift  term  can  be  expressed  by 


L = q S 


(105) 


where  Cj^  is  the  total  lift  coefficient,  q is  the  dynamic 


pressure,  and  S is  the  total  surface  area.  In  addition 


= + C,  a + C,  6e  + C.  ^ 

^o  S 2 V 


(106) 


where  C,  is  the  total  aircraft  lift  coefficient  for  a = 

^o 

6e  = 0,  Cy  is  the  total  aircraft  lift  curve  slope,  C.  is 

^6e 
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the  change  in  total  lift  coefficient  for  a unit  elevator  angle, 


C-  is  the  change  in  lift  coefficient  for  a change  in  pitch 


rate,  c is  the  mean  aerodynamic  chord  length,  and  is  the 


reference  pitch  rate.  For  the  problem  under  consideration. 


equals  zero.  Also,  since 


Qr  = ctj.  + Y 


then 


C,  = C,  + C.  6e_  + C, 

L L r Lr  r L 

a oe  q 


c(aj.  + y) 


2 V 


(107) 

(108) 


The  normal  acceleration  can  be  expressed  in  terms  of  curva- 
ture by 

(109) 


an  = V^K 


Substituting  the  above  relationships  into  equation  (104)  and 
rearranging  terms  results  in  the  following 


qS  6e 


(mV^K  + mg  cos  y - T e) 


qS 


"Lg  c (Ur  + Y ) 


2 V 


(110) 


It  is  now  necessary  to  estimate  the  reference  angle  of  attack 


and  its  derivatives.  This  estimated  angle  of  attack,  a^,  will 


not  be  used  explicitly  in  the  reference  command  generator,  but 
it  is  necessary  in  order  to  obtain  the  derivative  of  the  angle 
of  attack  for  use  in  equation  (110).  This  estimate  of  the 
angle  of  attack  derivative  will  also  be  used  to  calculate 
the  reference  pitch  rate  as  indicated  in  the  next  section. 
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In  order  to  arrive  at  an  estimate,  it  will  be  assumed  that 


“e 

a 


then 


Letting 


m(V^K  + g cos  Y) 

(qSC,  + T ) 
a 


(5SCl  T^) 


= V^K  + g cos  Y 


the  estimated  angle  of  attack  becomes 


(111) 


(112) 


(113) 


(114) 


a = 

e 

c*  V. 

A dg 

- (115) 

And  the 

derivative 

of  the  estimated  angle  of 

attack  is 

& 

e 

C.  V , + C.  V. 

A dg  A dg 

(116) 

or 

a 

e 

(0^/0^)%  V 

(117) 

where 

- t^) 

m 

(118) 

and 

3 V V K + K V* 

(119) 

The  second  derivatives  of  equations  (117),  (118),  and  (119) 
are  required  for  use  in  the  forthcoming  moment  equation  and 
can  be  evaluated  as 


'^dg  " 


(120) 


with 


- 2C,C,(q  SC,  + T^)  cUqSC,^  TJ 


m 
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m 


(121) 


And  = 3V2K+3VVK+5VVK+V2K  (122) 

dg 


where  the  dynamic  pressure  is 


q = 4 P V- 


(123) 


and  its  derivatives  are 


q = p V V 


(124) 


q = p (V  V + V2) 


(125) 


where  air  density,  p,  is  assumed  constant.  The  estimated 


thrust  derivatives  and  will  be  approximated  by 


first  backward  difference  equations.  The  expression  for 
K is  shown  in  equation  (102).  The  second  derivative  is 

• • 0 

cos^Y  - 4h'''  h’  K cos''y 


V K 

K = + V^i 


V 

6 K K h' 


- 3K2h' ' cos  y) 


(126) 


For  the  cubic  spline  the  fourth  derivative,  h'''',  will  be 
approximated  by  zero. 

Using  the  basic  relationship  of  equation  (108),  there 


still  remains  the  two  unknowns  a and  fie  . This  requires 

r r 


an  examination  of  the  moment  equation  in  order  to  have  two 
equations  and  two  unknowns.  The  moment  equation  can  be  ex- 
pressed by 


^yy**  ~ “a 


X* 


(127) 


: I 
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where  lyy  is  the  moment  of  inertia  about  the  y-axis,  11^ 
is  the  moment  due  to  the  aerodynamic  forces,  and  M,j,  is  the 
moment  due  to  thrust.  M,p  for  the  aircraft  under  study  is 
zero.  Thus,  the  total  moment  is  equal  to  the  moment  due  to 
the  aerodynamic  forces.  Expanding  the  right  side  of  the 
equation  yields 


“a  ' =11  5 ® <= 


(128) 


where  Cj^  is  the  total  aircraft  pitching  moment  coefficient, 
This  coefficient  can  be  expressed  as 

c q,. 


M M r M r M 

a oe  q 


2V 


• (129) 


where  Cj^  is  the  total  aircraft  pitching  moment  coefficient 
“cx 

versus  angle  of  attack  slope,  C„  is  the  change  in  total 

“6e 

aircraft  pitching  moment  coefficient  for  a unit  elevator 


angle 

, and  C,,  is  the  change  in  the  moment 

q 

coefficient  for 

unit 

change  in  pitch  rate. 

Since 

* • 

a = a + V 

^r  e ’ 

(130) 

then 

q = a + Y 

^r  e ' 

(131) 

Substituting  the  above  equations  into  equation  (127)  and 
rearranging  terms,  the  moment  equation  becomes 


C„  a + C„  6e 
Mr  M . 1 

a Se 


I (a  + y) 
yy  e ' 

q S c 


2 V 


(132) 


In  order  to  simplify  the  required  manipulations,  the  fol- 
lowing quantities  are  defined 
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Cr  + 

qS 

(mV^K  + mg  cos  y - T^c)  “e 


(133) 


(134) 


1 ( o ■*■  y)  C,,  (a  * y)  c 

yy  e ' _ Mq^^e  '' 

q S c 2 V 


(135) 


The  lift  equation  (equation  (110))  and  the  moment  equation 
(equation( 132) ) can  be  expressed  in  matrix  form  as 


a 6e 


The  leading  matrix  is  non-singular;  thus 


S ■ Si  -C 

a 6e  a 6e  M 


Solving  for  a^.  and  6e^  yields 


c„  ac.  — C. 

^6e  ^d  ^6e  “d 


'"“se  ■ 


g d g d 


" S ^L. 

g 6e  g 6e 


■ (136) 


g d 


(137) 


(138) 


(139) 


Now  that  the  reference  angle  of  attack  has  been  established, 
it  is  posible  to  update  the  estimated  first  derivative  of  the 
angle  of  attack,  g^,  for  use  in  defining  the  pitch  rate  ref- 
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erence,  q^. 


Reference  Pitch  Rate.  In  order  to  increase  the  accuracy 
of  the  reference  pitch  rate, a first  order  Taylor  series  ex- 
pansion is  accomplished  for  the  estimated  first  derivative 
of  the  angle  of  attack.  This  updated  derivative  is  then 
summed  with  the  first  derivative  of  the  flight  path  angle  to 
form  the  reference  pitch  rate. 

A first  order  Taylor  series  expansion  of  a results  in 
• • 


r e 


Aa 


3a 


(140) 


Using  Og  as  defined  in  equation  (117),  the  expansion  becomes 


r e 


CA(q  S + T^)  (Oj.  - Og) 

q 


m 


(141) 


The  updated  reference  pitch  rate  is 


Qj.  = Oj.  + Y 


(142) 


Reference  Thrust  Command.  The  reference  thrust  command 
can  be  calculated  from  the  drag  equation. 


mV  = T^  cos  (a^  + e)  - D - mg  sin  y 

X T ^ t 


(143) 


where  D is  the  total  drag  and  T^.  is  the  reference  thrust, 
From  equation  (92) 


V = - g sin  Y 


(144) 


Thus,  the  reference  thrust  is 

D 


cos(aj,  + e) 


(145) 
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Drag  can  be  expressed  by 


D = q S Cjj 


(146) 


The  total  drag  coefficient,  C^,  can  be  approximated  by 

r%  2 


= c, 


D. 


n A e 


(147) 


where  C.  is  the  total  aircraft  lift  coefficient,  is  the 

°o 

aircraft  drag  coefficient  at  Cj^  = 0,  A is  the  aspect  ratio, 
and  e is  the  efficiency  factor.  Since  the  reference  pitch 
rate  has  been  determined,  the  coefficient  of  lift  can  be 
evaluated  using  equation  (106).  Thus,  the  reference  thrust 
command  can  be  obtained  through  the  solution  of  equations 
(145),  (146),  and  (147). 


Reference  State  and  Control  Vectors . Utilizing  the  values 
obtained  thus  far  for  the  reference  angle  of  attack,  the  first 
derivative  of  the  angle  of  attack,  the  reference  thrust  com- 
mand, the  first  derivative  of  the  flight  path  angle,  and 
the  flight  path  angle,  the  reference  state  vector  becomes 


- - 

— 

' u 

V . cos  a 

r 

d r 

V,  sin  a 

r 

d r 

• • 

'Ir 

“r  Yd 

S 

9r 

“r  ^ Yd 

fie,. 

r 

r 

T 

T 

r 

r 

h 

r 

d 

(148) 


The  contol  vector  is  composed  of  the  elevator  and  thrust 


reference  states 


(149) 


The  proposed  system  can  now  be  implemented  as  described  in 
Chapter  I. 

Hybrid  Simulation  of  the  Terrain  Following  Task 

The  terrain  following  task  was  simulated  on  a Xerox 
Sigma  7,  32  bit,  112K  hybrid  computer  with  15  bit  plus  sign 
analog  to  digital  and  digital  to  analog  converters.  The 
command  follower  system  as  shown  in  Figure  1 of  Chapter  I 
was  implemented  on  the  hybrid  computer  using  the  reference 
command  generator  developed  previously  in  this  chapter,  the 
non-linear  aircraft  model  developed  in  Appendix  A,  and  the 
discrete  optimal  feedback  gain  matrix,  K,  based  on  the  per- 
formance index  developed  in  Chapter  III.  A short  segment  of 
terrain  (50,000  feet)  was  used  in  the  simulation  from  which 
a 100  foot  clearance  path  and  optimal  cubic  spline  path  were 
generated  (Ref  11).  The  sample  terrain  and  associated  paths 
are  shown  in  Figure  30.  The  desired  altitude  and  slope  at 
1000  foot  discrete  horizontal  range  increments  were  stored 
for  optimal  path  reconstruction  by  the  REFCOM  program.  The 
desired  path  was  constrained  to  remain  within  a minus  one  to 
a plus  two  "g"  normal  acceleration  boundary.  The  system  was 
observed  at  sample  rates  from  100  hertz  down  to  2 hertz  (twice 
the  Nyquist  rate).  The  state  response  of  the  system  for  100, 
5,  and  2 hertz  sample  rates  is  shown  in  Figures  31  through  38. 
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Figure  30.  Optimal  cubic  spline  path  and  terrain  for  hybrid  simulation. 


The  error  vector  quantities  (state  minus  reference  vector) 
are  shown  in  Figures  39  through  45. 

As  shown  in  the  associated  figures,  the  state  response 
is  nearly  identical  for  sample  rates  from  100  to  5 hertz.  The 
primary  difference  is  in  the  response  of  the  elevator  state 
as  indicated  in  Figure  35.  A greater  maximum  and  minimxim 
deflection  is  required  at  the  five  hertz  sample  rate.  It 
is  interesting  to  note  the  state  behavior  at  a sample  rate 
of  two  hertz.  Only  slight  degradation  of  the  state  response 
can  be  observed.  In  fact,  the  elevator  response  has  been 
smoothed  at  this  low  sample  rate. 

Of  further  interest  is  the  relationship  shown  when  the 
thrust  time  history,  Figure  36,  is  compared  to  the  angle  of 
attack  history.  Figure  37.  Since  the  reference  command  gen- 
erator was  based  on  a simplified  constant  energy  control,  as 
angle  of  attack  is  increased  depleting  the  energy  level 
through  aircraft  "g"  loading  a resulting  thrust  command  is 
required  to  restore  the  depleted  energy.  To  further  reduce 
thrust  fluctuation,  in  order  to  achieve  less  engine  cycling 
and  wear,  it  would  be  necessary  to  incorporate  a means  of 
allowing  the  energy  level  to  vary  between  a lower  and  upper 
limit.  This  could  be  implemented  through  software  changes 
quite  readily. 

Even  though  the  altitude  response  as  shown  in  Figure  38 
is  remarkably  similar  for  a wide  range  of  sample  rates,  it 
is  necessary  to  examine  the  error  between  actual  and  reference 
altitude  to  determine  acceptability.  This  error  is  illus- 
trated in  Figure  39.  The  maximum  absolute  error  is  13  feet 
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at  100  hertz  and  increases  to  22  feet  at  a five  hertz  sample 
rate.  For  a two  hertz  rate,  the  absolute  error  increases  to 
34  feet.  Thus,  the  five  hertz  sample  rate  could  be  deemed  to 
be  quite  acceptable  in  this  tracking  task,  while  the  two  hertz 
system  could  be  deemed  as  aedequate.  This  is  also  verified 
by  the  other  state  variable  responses  at  the  two  hertz  sample 
rate.  If  it  were  required  to  guarantee  a minimum  clearance 
path  altitude,  the  optimal  path  altitude  could  be  shifted  up- 
ward by  the  maximum  error  to  insure  achievement  of  the  desired 
clearance  minimum.  It  is  also  entirely  possible  to  decrease 
the  error  quantity  by  tuning  of  the  controller. 

The  other  error  states  are  included  in  Figure  40  through 
45  for  the  interested  reader. 
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31.  Longitudinal  velocity 


as  a function  of  time 


2.  Vertical  velocity,  w,  as  a function  of  time 
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34.  Pitch  angle 


as  a function  of  time 


Figure  35.  Elevator  angle,  6e,  as  a function  of  time 


Figure  36.  Thrust,  T,  as  a function  of  time 


Figure  37.  Angle  of  attack,  a,  as  a function  of  time 
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Altitude,  h,  as  a function  of  time 
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Figure  40.  Error  in  longitudinal  velocity,  Au,  between 
actual  and  reference. 


Figure  41.  Error  in  vertical  velocity,  Aw,  between  actual 
and  reference. 
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Figure  42.  Error  in  pitch  angle,  A0,  between  actual  and 
reference. 


Figure  44.  Error  in  elevator  deflection  angle,  A6e,  between 
actual  and  reference. 
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Figure  45.  Error  in  thrust  setting,  aT,  between  actual 
and  reference. 
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VI . Conclusions  and  Recommendat ions 

Conclusions 

The  breakthroughs  oi  the  recent  past  have  provided  the 
control  engineer  with  a new  method  for  the  implementation  of 
an  aircraft  flight  control  system.  Through  the  use  of  state 
space,  continuous  optimal  control,  and  discrete  optimal  con- 
trol, a digital  flight  control  system  has  been  designed  for 
the  terrain  following  task.  After  formulating  the  aircraft 
linear  perturbation  model,  the  deterministic  regulator  prob- 
lem was  solved  iteratively  to  provide  the  state  cost  and 
control  cost  matrices  for  the  desired  continuous  closed  loop 
system.  The  selection  of  these  matrices  was  based  on  the 
Cornell  Aeronautical  Laboratory  "thumbprint"  which  specified 
the  short  period  natural  frequency  and  damping  ratio.  Once 
the  matrices  were  selected,  an  analog  computer  simulation 
was  accomplished  to  verify  the  "goodness"  of  the  selection. 

As  a side  note,  the  controller  was  modified  to  yield  a C* 
response  for  comparison  with  C*  criteria.  With  the  continu- 
ous cost  matrices  defined  for  the  desired  closed  loop  root 
location,  the  system  was  discretized,  as  was  the  continuous 
quadratic  performance  index,  to  form  a discrete  deterministic 
regulator  problem. 

The  discrete  regulator  problem  was  solved  as  a function 
of  sample  rate  to  determine  a minimum  rate  for  sampling. 

This  direct  digital  design  technique  has  the  advantage  of 
not  introducing  errors  into  a discrete  system  analysis  since 
the  zero  order  hold  is  accounted  for  precisely.  The  effects 
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of  sample  rate  on  the  system  were  then  examined.  The  finite 
word  length  of  a digital  system  was  shown  to  produce  the 
equivalent  of  a noise  being  introduced  into  the  system.  A 
sample  rate  of  five  hertz  was  shown  to  be  high  enough  to 
aedequately  form  the  desired  controls.  A reference  command 
generator  was  developed  using  as  a basis  a form  of  constant 
energy  control.  A hybrid  simulation  was  then  completed  using 
a reference  path  generated  by  Funk's  optimal  cubic  spline 
algorithm.  The  aircraft  was  shown  to  track  the  desired  path 
in  a highly  acceptable  manner. 

Recommendations  for  Further  Study 

Since  the  system  studied  was  assumed  to  be  a perfect 
system  with  body  bending  and  sensor  noise  not  considered,  it 
is  recommended  u.  stochastic  study  of  the  problem  be  conducted  ' 

with  these  items  taken  into  consideration.  Also  of  interest  ' 

would  be  the  incorporation  of  an  estimator  or  observer  since  i 

for  this  study  all  states  were  assumed  available.  The  noise 

i 

due  to  finite  word  length  could  also  be  taken  into  account 

using  a stochastic  approach.  | 

Another  realm  of  interest  would  be  an  assessment  of  the  ! 

I 

required  memory  size  for  the  control  system.  This  study  I 

j 

should  include  the  requirements  for  the  optimal  cubic  spline  ' 

algorithm. 

As  a final  recommendation,  the  approach  taken  in  this 
study  could  be  adopted  for  tasks  other  than  terrain  following. 

For  instance,  weapon  delivery  or  air  to  air  tracking  would 
be  of  great  interest. 

82  I 

J 


B^bliographj 


Hartmann,  G.  L.  and  C.  A.  Harvey.  "Digital  Adaptive  Con- 
trol Laws  for  the  F-8."  American  Institute  of  Aeronautics 
and  Astronautics.  76-1952;  319-330  (August  1976). 

A-7D  Digital  Multimode  Flight  Control  System  Flight  Test 
and  Weapon  Delivery  Evaluation . AFFDL-TR-76-121 . Wright- 

Patterson  Air  Force  Base,  Ohio:  Air  Force  Flight  Dynamics 
Laboratory,  December  1976. 

Advanced  Fighter  Digital  Flight  Control  System  (DFCS)  Def- 
inition Study.  W0728-FR.  Minneapolis,  Minnesota:  Honey- 
well, Inc.,  March  1975. 

Powell,  J.  D. , E.  Parsons,  and  M.  G.  Tashker.  Aircraft 
Digital  Control  Design  Methods . SUDAAR  500.  Stanford, 
California;  Stanford  Uniyersity  Guidance  and  Control 
Laboratory,  Department  of  Aeronautics  and  Astronautics, 
Feruary  1976. 

Katz,  P.  and  J.  D.  Powell.  Selection  of  Sampling  Rate 
for  Digital  Control  of  Aircrafts . SUDAAR  486.  Stanford, 
California:  Stanford  Uniyersity  Guidance  and  Control 
Laboratory,  Department  of  Aeronautics  and  Astronautics, 
September  1974. 

Tobie,  H.  N. , E.  M.  Elliott,  and  L.  G.  Malcom.  A New 
Longitudinal  Handling  Qualities  Criterion . Seattle, 
Washington:  Commercial  Airplane  Division,  The  Boeing 
Company,  1965. 

Roskam,  J.  Flight  Dynamics  of  Rigid  and  Elastic  Airplanes , 
Part  1.  Kansas:  Roskam  Ayiation  and  Engineering  Corpora- 
tion, 1976. 

Funk,  J.  E.  "Optimal-Path  Precision  Terrain  Following 
System."  American  Institute  of  Aeronautics  and  Astro- 
nautics , 76-1958:  383-390  (August  1976). 

Bryson,  A.  E.  and  W.  E.  Hall.  Optimal  Control  and  Filter 
Synthesis  by  Eigenyector  Decomposition . SUDAAR  436.  Stan- 
ford, California:  Stanford  Uniyersity  Department  of 
Aeronautics  and  Astronautics,  Noyember  1971. 

Kisslinger,  R.  L.  and  M.  L.  Wendl.  Suryjyable  Flight  (Con- 
trol System  Interim  Report  No . ^ Studies , Analysis  anH 
Approach . AFFDL-TR-71-20  Supplement  1.  Wright-Patterson 
Air  Force  Base,  Ohio:  Air  Force  Flight  Dynamics  Laboratory, 
May  1971. 

Funk,  J.  E.  Terrain  Following  Control  Based  on  an  Optimized 
Spline  Model  of  Aircraft  Motion . ASD/XR-TR-75-22 . Wright- 


83 


f 


Bibliography  (Continued) 


Patterson  Air  Force  Base,  Ohio:  Deputy  for  Development 
Planning,  Aeronautical  Systems  Division,  November  1975. 


I 


84 


r 

Appendix  A 

f Non-Linear  Aircraft  Model 

For  the  problem  under  study,  it  is  assumed  that 

manue- 

vering  is  restricted  to  the  longitudinal  plane;  thus 

» the 

longitudinal  equations  of  motion  become 

m(U  + WQ)  = - mg  sin  0 + F + F 

^x  ^x 

(A-1) 

• 

m(W  - UQ)  = mg  cos  9 + Fg 

z 

(A-2) 

lyyQ  = “a  ^ “t 

(A-3) 

Further,  for  wings  level  non-turning  flight 

Q = 6 

(A-4) 

The  approximations  for  the  elevator  and  throttle  actuators 

■* 

are 

6e  = - 10  6e  + 10  6e^ 

(A-5) 

• 

and  T = - T + 

c 

(A-6) 

Altitude,  h,  can  be  defined  by 

h = V sin  Y 

(A-7) 

The  force  produced  by  thrust  in  the  x-direction  is 

F,j,  = T cos  e 

(A-8) 

where  e is  the  angle  the  engine  thrust  line  makes  with  the 

x-axis.  In  a similar  fashion 

F,j,  = - T sin  c 

z 

85 

(A-9) 

J 

For  the  aircraft  under  study,  e is  approximately  2.9  degrees. 
The  thrust  vector  acts  through  the  center  of  gravity;  thus, 
the  moment  due  to  thrust  is  zero. 

The  aerodynamic  forces  acting  in  the  x-direction  can  be 
expressed  by 

F = - D cos  a + L sin  a (A-10) 

^x 

the  total  drag  is  represented  by 

D = Cp  q S (A-11) 

where  is  the  total  aircraft  drag  coefficient,  S is  the 
surface  area,  and  q is  the  dynamic  pressure  acting  on  the 
aircraft.  Dynamic  pressure  can  be  evaluated  by 

q = 4 p V*  (A-12) 

where  p is  the  atmospheric  density  and  V is  the  total  vel- 
ocity. The  total  velocity  is  defined  by 

= U*  + W*  (A-13) 


The  atmospheric  density,  p,  is  assumed  constant  over  the 
realm  of  flight  and  is  evaluated  at  sea  level  ( P^j^  = .002377 
slugs/ft ^).  The  total  aircraft  drag  coefficient  will  be 
approximated  by 


C 


D 


•ff  A e 


(A-14) 


where  C.  is  the  total  aircraft  lift  coefficient,  is  the 
aircraft  drag  coefficient  at  Cj^  = 0,  A is  the  aspect  ratio, 
and  e is  the  efficiency  factor.  It  is  further  assumed  that 
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the  Mach  and  Reynolds  numbers  are  approximately  constant 
for  the  flight;  thus,  the  total  aircraft  drag  coefficient 
can  be  approximated  at  this  Mach  nvunber. 

The  aerodynamic  forces  in  the  z-direction  can  be  summed 

into 

F = - L cosa  - D sina  (A-15) 

Lift  can  be  expressed  by 


L = q S 


(A-16) 


The  lift  coefficient  will  be  approximated  by 


c Q 

Cj^  — + Cj^  a + Cj^  6e  + Cj^  — 

o a 6e  q 2 V 


(A-17) 


Where  is  the  total  aircraft  lift  coefficient  for  a=  6e=  0 

(for  this  aircraft  C. 4=  0),  is  the  total  aircraft  lift- 

■L’—  L 

o a 

curve  slope,  C,  is  the  change  in  total  aircraft  lift  coeffi- 
^6e 

cient  for  a unit  elevator  angle,  C.  is  the  change  in  lift 
coefficient  for  a change  in  pitch  rate,  and  c is  the  mean 
aerodynamic  chord  length. 

The  moment  due  to  the  aerodynamic  forces  is 


M 


A 


Cm  q S a 


(A-18) 


where  Cy  is  the  total  aircraft  pitching  moment  coefficient. 
This  coefficient  will  be  approximated  by 


o a 6e 


iL£. 


a 6e  "q  2 V 

where  is  equal  to  zero  for  this  aircraft,  C 


M 


(A-19) 
is  the 


o a 

total  aircraft  pitching  moment  coefficient  versus  angle  of 
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r- 

attack  slope,  C,,  is  the  change  in  total  aircraft  pitching 
“6e 

moment  coefficient  for  a unit  elevator  angle,  and  is  the 

q 

change  in  the  moment  coefficient  for  unit  change  in  pitch 
rate. 


The  coefficient  values  and  initial  conditions  for  the 
non-linear  simulation  are 


f 

f 

= 

2.066 

c = 

18.6  ft 

a 

— 

.1383 

A 

2.5 

= 

-.542 

e = 

.9 

q 

"-0 

Cm 

s 

.0106 

-.1833 

'yy  = 
s 

.38  X 10®  slug  ft 
647  ft^ 

a 

Cm 

6e 

s 

-.1529 

p 

.002377  slugs/ft® 

Cm 

* 

-.651 

g 

32.174  ft/sec*^ 

"q 

^o 

= 

42.65  ft/sec 

e = 

2.9  degrees 

k 

Vo 

s 

882.0  ft/sec 

m = 

1709.4  slugs 

Qo 

0.0  deg /sec 

- 3.324  degrees 

®o 

= 

2.772  degrees 

“o  = 

2.772  degrees 

«o 

= 

0.0  feet 

To  = 

7091.2  pounds 

«o 

= 

300  feet 

Uo 

= 

880.96  ft/sec 

J 
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The  following  program  requires  precomputed  storage  of 
the  desired  cubic  spline  path  altitude  and  slope  at  discrete 
range  intervals.  Subroutine  PATHGEN  reconstructs  the  required 
path  and  forms  its  derivatives  for  a given  range.  The  re- 
quired path  slope  and  altitude  was  provided  by  the  optimized 
spline  algorithm  developed  by  Funk  (Ref  11).  Certain  other 
parameters  are  required  to  interactively  operate  the  program 
from  an  intercom  terminal.  The  complete  digital  program  for 
the  hybrid  simulation  was  not  included  due  to  its  length. 
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Southeast  Asian  tour,  he  completed  a four  year  tour  at  RAF 
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